Heat capacity estimators for random series path-integral methods by finite-difference 

schemes 
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Previous heat capacity estimators used in path integral simulations either have large variances that 
grow to infinity with the number of path variables or require the evaluation of first and second order 
derivatives of the potential. In the present paper, we show that the evaluation of the total energy by 
the T-method estimator and of the heat capacity by the TT-method estimator can be implemented 
by a finite difference scheme in a stable fashion. As such, the variances of the resulting estimators 
are finite and the evaluation of the estimators requires the potential function only. By comparison 
with the task of computing the partition function, the evaluation of the estimators requires k + 1 
times more calls to the potential, where k is the order of the difference scheme employed. Quantum 
Monte Carlo simulations for the Nei3 cluster demonstrate that a second order central-difference 
scheme should suffice for most applications. 
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I. INTRODUCTION 

It is said that path integral methods transform a quan- 
tum equilibrium problem into a classical one by judi- 
cious use of dimensionality^ Yet, the computation of 
the average energy 2 : 3 : 4 : 5 : 6 : 7 : 8 : 9 : 10 : 11 : 12 : 13 : 14 or the heat 
capacitjii£*i£ of a quantum canonical ensemble reveals 
that the quantum-classical analogy is far from being triv- 
ial, even if distinguishable particles are assumed. One 
observes an increase in the computational time not only 
with the number of path variables considered, but also 
with the dimensionality of the system. This is so because 
estimators of finite variance usually involve first or sec- 
ond order derivatives of the potential. The number of 
such derivatives scales linearly or quadratically with the 
number of degrees of freedom of the system. For exam- 
ple, numerical studies of even moderately large quantum 
clusters are severely hindered by this substantial increase 
in the number of quantities that must be evaluated^ 

Recently, Predescu and DolU4 have observed that 
a simple rescaling of the Brownian bridge entering 
the Feynman-Kac formulatf^S*^ produces path-integral 
techniques for which the dependence with the tempera- 
ture of the path distributions is buried into the potential 
part of the imaginary-time action. Formal differentia- 
tion of the logarithm of the partition function leads to a 
special form of the thermodynamic estimator (T-method 
estimator) that does not have the variance difficulties as- 
sociated with the Barker estimator for large numbers of 
path variables Even though the resulting T-method 
estimator closely resembles the virial estimators j&iSii 3 . it 
does not rely on the virial theorem to recover the kinetic 
energy. For instance, this T-method estimator produces 
correct results even for potentials that are not confined, 
e.g. a free particle. Therefore, the variance of the T- 



method estimator is lower than that of the virial estima- 
tor because the classical part of the energy is explicitly 
introduced as a constant and is not obtained from the 
virial theorem. In a recent study of the (112)22 cluster 
at the temperature of 6 K, 20 difficulties associated with 
the virial estimator for low-temperature systems 4 * 7 ^ were 
not observed to appear for the T-method estimator in- 
troduced by Predescu and Doll. Such differences between 
the estimators are even more significant for the heat ca- 
pacity estimators and will be revealed in the present pa- 
per by comparing the statistical errors for the Predescu 
and Doll-type estimators with those for the double virial 
estimator 

In order to avoid any confusion with earlier estima- 
tors, we mention that in the present article by T-method 
and H-method estimators we understand the respective 
ener gy e stimators introduced by Predescu and Doll in 
Ref. Hi By TT-method and TH-method estimators, 
we understand the heat capacity estimators that are ob- 
tained from the corresponding energy estimators by tem- 
perature differentiation. 

The T-method estimator is closely related and simi- 
lar in form to the centroid virial estimator There 
are however two differences. First, the T-method estima- 
tor involves fluctuations of the Brownian bridge relative 
to one arbitrary point. The centroid virial estimator in- 
volves similar fluctuations but relative to the path cen- 
troid. It can be shown that the ratio between the aver- 
age square fluctuations of the Brownian bridge relative 
to some preferential point and to the path centroid is 
2^1 Thus, the two estimators have similar behavior with 
the nature of the quantum system, the temperature, and 
the Monte Carlo sampling method, though the centroid 
virial estimator may exhibit a slightly lower variance. 

A second and more important difference, which con- 
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stitutes the starting point of the present development, is 
the fact that the T-method estimator is a veritable ther- 
modynamic estimator, in the sense that it is obtained 
by temperature differentiation of the quantum partition 
function (however, as discussed in a previous paragraph, 
one needs to utilize a special form for the Feynman-Kac 
formula, with the temperature dependence of the path 
distribution buried into the potential). The tempera- 
ture differentiation can be implemented numerically by 
a finite-difference scheme and leads to numerically stable 
algorithms that do not require derivatives of the poten- 
tial. This observation proves to be extremely important 
for heat capacity calculations because formal tempera- 
ture differentiation leads to expressions involving all first 
and second order derivatives of the potential. By numeri- 
cal temperature differentiation, one obtains an important 
speed-up in the evaluation of the above mentioned ther- 
modynamic properties, especially for large dimensional 
systems or for complicated potentials. 

In this article, we also propose an analytical heat ca- 
pacity estimator that involves the first derivatives of 
the potential only. This is obtained from the analytical 
form of the TT-method estimator by an integration by 
parts suggested by Predescu and Doll in the derivation 
of their special H- method energy estimator iH The two 
estimators, called in this paper the TT-method estima- 
tor and the modified TT-method estimator respectively, 
may have slightly different variances. As discussed in the 
previous paragraph, the first one is to be implemented by 
finite-difference, whereas for the second one we shall use 
exact analytical formulas. 

The relative merits of the new heat capacity estima- 
tors will be demonstrated for the Nei3 cluster. For this 
example, we provide a comparison of the statistical er- 
rors of the new estimators with those of the double virial 
estimator utilized by Neirotti, Freeman, and Dolli^ We 
shall also clarify a number of issues raised in the Neirotti, 
Freeman, and Doll study of this neon cluster. The nu- 
merical simulation presented serves to demonstrate the 
power of the path integral approach utilized as well as 
to provide essentially exact numerical data necessary for 
comparison in the development of quantum approxima- 
tions that can be employed for larger or more complicated 
systems 



II. THERMODYNAMIC ENERGY AND HEAT 
CAPACITY ESTIMATORS 

In this section, we derive the formal expressions for 
the heat capacity of a d-dimensional canonical quantum 
mechanical system made up of distinguishable particles. 
The particles have masses {mo,i; 1 < i < d} and move in 
the potential V(x). The vector x, the transpose of which 
is x T = (xi, . . . ,Xd), denotes the position of the parti- 
cles in the configuration space M. d . The canonical system 
is characterized by inverse temperature ft = l/(ksT). 
Its average energy and heat capacity can be obtained 



by temperature differentiation of the partition function 
Z{ft), producing the formulas 



{E)\ 



I dZ(/3) 
Z(ft) d/3 



and 



{Cv) T p T = k B 



ft 2 d 2 Z(j3) 
Z(f3) d/3 2 



ft dZ(ft) 
Z(ft) dp 



(1) 



(2) 



respectively. The partition function of the system is ob- 
tained as the integral over the configuration space of the 
diagonal density matrix 



Z(ft) 



p(x; ft)dx.. 



(3) 



In the path-integral approach, the density matrix is 
evaluated with the help of the Feynman-Kac formula. 
We split the present section into two parts. In the first 
part, we shall discuss the random series implementation 
of the Feynman-Kac formula and introduce some relevant 
notation. In the second part, we deduce the formal ex- 
pression of the TT-method heat capacity estimator and 
discuss its numerical implementation by finite-difference 
schemes. Then, we derive the modified TT-method esti- 
mator, the analytical expression of which involves first- 
order derivatives of the potential only. 



A. Random series path integral techniques 

In the random series implementation of the Feynman- 
Kac formula, the density matrix of a one-dimensional 
quantum system is obtained as follows. 14 Let {Afe(T)}fe>i 
be a system of functions on the interval [0, 1] that, to- 
gether with the constant function Ao(t) = I, make up an 
orthonormal basis in L 2 [Q, I]. Define 



A fe (i) 



\k(u)du. 



Let n denote the space of infinite sequences a 
(ax, &2, ■ • •) and let 



dP[d] = [] Mak) 



(4) 



fc=i 



be the probability measure on O such that the coordi- 
nate maps a —> afc are independent identically distributed 
(i.i.d.) variables with distribution probability 

dn{ai) = ^Le~ a - /2 da t . (5) 



/2tt 

Then, the Feynman-Kac formula reads 



p(x,x';ft) 
p fp (x,x';ft) 



dP[a]exp<^ -ft / V x r (u) 



du 



k=l 



(6) 
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where x r (u) = x + (x' — x)u and a = (h 2 (3/mo) 1 ^ 2 . The 
quantity pf p (x,x';(3) represents the density matrix for a 
similar free particle. The series 



B° u (a) 



k=l 



a fe A fc (u) 



represents a stochastic process equal in distribution to a 
standard Brownian bridge. 

For a d-dimensional system, the Feynman-Kac for- 
mula is obtained by employing an independent random 
series for each additional degree of freedom. As such, 
we consider the space £l d made up of all sequences 
a = (ai, a 2 , . . .) of vectors 



a/, 



a>d,k 



and denote the line i of a by aj — (ffl^o, a i,ii ■ ■ ■)■ O n the 
space il d , we define the probability measure 



dP[a] =Y[dP[oi] (7) 



with 



dP[ai] = JJ d/i(oi,fc). 
fc=i 

Under this probability measure, the coordinate maps 
a ~~ * Oi,jfe a re i-i.d. standard normal variables. We also 
consider the vector a T = (<7i , . . . , a^) of components 
(Ti = (K 2 p/moj) 1 ' 2 and let x r (u) = x + (x' — x)u be 
a straight line connecting the points x and x'. Then, the 
Feynman-Kac formula reads 



p(x,x';/3) 



dP[a]exp^ - f3 / V x r (u) 



cr^a fc A fc (w) 



fc=i 



(8) 



where 



cra fc 



&d0.d,k 



The series 



^a fc A fe (u) 
fc=i 



is equal in distribution to a d-dimensional Brown- 
ian bridge (a vector- valued stochastic process whose 
components are independent one-dimensional Brownian 
bridges) . 



In practical applications, one cannot work with the full 
random series implementation of the Feynman-Kac for- 
mula. Instead, one considers finite-dimensional approxi- 
mations to Eq. © , the simplest of which have the general 
form 



p n (x,x';/3) 
P/ P (x,x';/3) 



dP %] 



<\\]> ) - P I V 



qn+p 



x r (u)+a ^ a fc A„ !fc (u) du>, (9) 



fc=l 



where q and p are some fixed integers. The functions A„ t k 
are chosen so that to maximize the rate of convergence 

p n (x,x';/3) -> p(x,x';/3). 

Though the problem of maximizing the order of conver- 
gence is still far from a final resolution, several schemes 
in the larger class of reweighted techniques were proven 
to have cubic or quartic asymptotic convergence^^ 
The construction of the functions A rit k and of associ- 
ated quadrature schemes for the computation of the 
path averages appearing in Eq. have been dis- 
cussed elsewhere For the numerical examples 
presented in this article, we use a so-called Levy- 
Ciesielski reweighted path integral method having quar- 
tic convergence^ To a large extent, the analytical ex- 
pressions of the functions K nt k{u) and the nature of the 
quadrature schemes are not important for the present de- 
velopment. For more information, the reader is advised 
to consult the cited references. 

To simplify notation, we introduce several auxiliary 
quantities B° n (a), ?7 n (x, x', a; and X n (x, x', a; /3), 
defined by the expressions 



qn+p 



#°,n( a ) = Z2 a fcAn,fc(u), 



fc=l 



Z7„(x,x',a; 0) 



V 



x r (w) + aB° (a) 



du, 



(10) 



fill 



and 



X n (x, x', a; j3) = pf p (x, x'; (3) exp [-/?[/„ (x, x', a; /?)] , 

(12) 

respectively. The similar relations for the full Feynman- 
Kac formula are denoted by i?°(a), £/oo(x, x', a; f3), and 
Aoo(x, x', a; f3), respectively. With the new notation, 
Eq. (0 becomes 

p n (x,x';/3)= / dP[a]A„(x,x',a;/3), (13) 

whereas the Feynman-Kac formula reads 

p(x,x';/3)=/ dP[a]X 00 (x,x',a;/?). (14) 
Jn d 
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In this paper, we make the convention that x' is 
dropped whenever x = x'. In order to arrive at the defi- 
nition of the energy and the heat capacity estimators, it 
is convenient to introduce the quantities 

K n (x,a,/J,e)-—— — — — - e 
A„(x,a; fj) 

x exp [-(3eU n (K, a; pe) + PU n (x, a; p)] (15) 



and 

x exp [-^[/^(x, a; /3e) + PU^x, a; /?)] 
respectively. We have 



(16) 



and 



/3 dZ{(3) I R <idxJ nd dP[a}X Q0 {x 1 a;P)f t R 00 {x,a;P 1 e) 



Z{P) dp 



J Rd dx f nd dF[a]Xoo (x, a; P) 



(17) 



P 2 d 2 Z(P) lRddxJ nd dP[a\X 00 (x 1 a.;P)-^R 00 (x,a- 1 P,e) 



Z{P) dp 2 



Jr-* rfx /n<J d-P[a]^oo(x, a; P) 
I 



(18) 



The quantities above can be evaluated by Monte Carlo integration as the limit n — > oo of the sequences 
I 



and 



P dZ n (p) J Rd dxJ nd dP[a}X n {x,a;P)£R n (x,a;P,e] 



Z n {P) dP 



Jr<* dx In* dP[a]X n (x, a; P) 



(19) 



P 2 d 2 Z n (P) _ J Rcl dxf nd dP[a}X n (x,a;p)£rR n (x,a;P,e) 



z n (p) dp 2 



J Rd dx J Qd dP[a]X n (x, a; P) 
I 



(20) 



respectively. 

In the finite-difference scheme that is advocated in this 
paper, the derivatives against e appearing in Eqs. i|19|) 
and (|2U[) may be evaluated numerically by finite differ- 
ence. Such an approach is expected to be much faster 
than the analytical evaluation of the derivatives, es- 
pecially for large dimensional systems or systems with 
complicated potentials. Though higher order central- 
difference schemes can be employed, a second order 
scheme produces 



— R n (x,a;P,e) 
x [R n (x,a;P, 1 



feo) 



- Rn(x,a;P,l- e ) 



and 



d 2 

-j-^R n (x, a;/3,e) 



e 2 [Rn(x,a;P, 1 + e ) 



2i?„(x, a; P, 1) + R n (x, a; 0, 1 - e )] 



with error of order O(cq). Such a direct approach may 
prove useful for highly quantum systems or for patholog- 
ical systems, as for instance a particle in a box. However, 
for smooth enough potentials, the alternatives that are 
analyzed in the following subsection may prove to be su- 
perior. 



B. Expressions of the heat capacity estimators 

In this subsection, we shall put the relevant quantities 
entering the expression of the heat capacity estimator in 
a form that is exact in the high-temperature limit or in 
the limit that the physical system is classical. For this 
purpose, we assume that exp[— pV(x)] has second order 
Sobolev derivatives as a function of x. In the second 
part of the present subsection, we shall derive a special 
analytical expression for the heat capacity estimator that 
employs the first order derivatives of the potential, only. 



5 



This modified heat capacity estimator gives results iden- 
tical to the first one, but it may have a slightly different 
variance. 

By explicit computation, one argues that 



The expression given by Eq. (|24|l is not computation- 
ally very convenient because it requires the evaluation of 
d(d + l)/2 path averages for as many different second- 
order derivatives 



and 



de 



i?„(x,a;/3, e) 



~- 0U n (x, a; 0) - PjU n {^ a; 0e) 



(21) 



de 2 
d 



i?„(x, a;/3,e) 



de 



i?„(x, a;/3,e) 



2/3-C/„(x,a;/3e) 
2 de 



/3^2^n(x,a; 0e) 



:=1 



(22) 



The first and second derivatives of the function 
I/ n (x, a;/3e) around the point e = 1 can be evaluated by 
finite difference, as shown in the preceding subsection. 
However, we notice that in the limit that the physical 
system behaves classically, f7„(x, a;/3e) w V^(x) and the 
derivatives against e vanish. Moreover, in this limit any 
finite-difference scheme produces the exact classical re- 
sults. It is therefore apparent that the utilization of the 
derivatives of the functions U n (x, a;/3e) instead of the 
derivatives of i?„(x, a; fj, e) has certain numerical advan- 
tages, increasing the range of acceptable values for the 
discretization step eo- 

We now proceed and compute the analytical expression 
of the derivatives of the function f/„(x, a; (3e). We have 

d 



— Z7 n (x,a; (3e) 
de 



x + <„(a) i?°i(a)du, (23) 



It is for this reason that we advocate the use of a finite 
difference scheme instead of the analytical formulas. For 
large enough physical systems or for complicated poten- 
tials for which the derivatives are not readily available, 
the finite difference scheme will enjoy a significant com- 
putational advantage. Parenthetically, Eq. (|24|l shows 
that the TT-method heat capacity estimator is similar 
in form to the double virial heat capacity estimatosi^ or 
to the centroid double virial heat capacity estimator. 16 
However, it has the distinctive feature (characteristic of 
the Barker estimators) that it can be implemented by a 
finite-difference scheme, yet it maintains to a good degree 
the low variance of the centroid double virial estimator. 



where 



qn+p 
fe=l 



One also computes 

d 2 1 d 

-p[/„(x,a;/3e) ^ = —^<Ji 



x / diV 
Jo 



x + aBl n (a)] B° u %{H)du + - ^ a^j 



x / dijV 



x + <„(a) B^(a)B2^(a)du, (24) 



For strongly quantum systems, as for instance low tem- 
perature hydrogen or helium clusters, there is sometimes 
the need to validate the convergence of the path inte- 
gral methods by employing the agreement between the T- 
method and the H-method energy estimators. 20 As shown 
by Predescu and DolljA^ the H-method estimator can be 
put into the "force-force correlation" form by a simple 
integration by parts. This form requires the first order 
derivatives of the potential only. In such cases, given that 
the first order derivatives of the potential are computed 
anyway, it would be advantageous if we could evaluate 
the heat capacity as a functional of these derivatives only. 
This can actually be done (again by integration by parts) 
as follows. Observe that 
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Jo L J 

dx 3 e-^° v \* + ° B °^ du {3 ( [\v\k + *B°, n (&j\ I%\(&)B°* n (&)du \ { I 'djV [x + a<„(a)l du \ . (25) 



r 



Therefore, Eq. I|24ll can be replaced for the purpose of 
evaluating the heat capacity by 



de 2 



Un 



1 d 

(x,a;/3e) e--^ 

e— 1 4 ' — rf 

i=l 



x / a 4 v 

Jo 



x + <rB° >fl (a) B2',n(a)d« + J 



x + oBl n (a)l (a)fl°£ (a)du |> ( 2( i ! 



We utilize the sign of equivalence = in the relation above 
to warn the reader that the equality implied by Eq. I|26|l 
does not hold in the strict sense. Rather, it means 
that the expression to the right of the sign = produces 
estimates identical to the ones obtained by employing 
Eq. J2U, though it may exhibit a different variance. 
The resulting heat capacity estimator will be called the 
modified TT-mcthod estimator and will be denoted by 
(Cy)™ TT henceforth. Eq. (|2l)|l still involves d 2 path aver- 
ages to be computed (which may become prohibitive for 
large dimensional systems), but this time the averages 
involve quantities that are computed anyway. Expensive 
calls to the second order derivatives of the potential are 
avoided. 



III. A NUMERICAL EXAMPLE 

We shall test the merits of the two heat capacity esti- 
mators discussed in the previous paragraph on a cluster 
of N p = 13 neon atoms using a special path integral tech- 
nique introduced in Ref. |2J and having quartic asymp- 
totic convergence with respect to the number of path 
variables. The numerical implementation of this method 
is similar to the Levy-Ciesielski reweighted method uti- 
lized in Ref. |2fJ and will not be reviewed here. 

Quantum studies of small Lennard-Jones neon clusters 
(N p < 7) by ground-stateSSi2LS& or finite-temperature 
methods22i2£ have revealed that the quantum effects are 
important, leading to large zero-point energies. By com- 
parison, studies of larger clusters are relatively scant. 
The Nei3 cluster is interesting because it is the smallest 
Lennard-Jones cluster that presents an effective classical 
melting point (at about 10 K) marking a transition from 



a rigid to a liquid-like phase. The pronounced quantum 
effects have been found to lower the transition tempera- 
ture by about 1 Ki 15 i 31 However, quantum heat capaci- 
ties reported in literature and computed by path integral 
methodsi^i or semiclassical techniques^ are not suf- 
ficiently accurate due to large statistical or systematic 
errors. To demonstrate the advantage of the new esti- 
mators, we propose to compute the heat capacity of the 
Nei3 cluster over the range of temperatures 4— 14 K, with 
a statistical error (defined in the present article as two 
times the standard deviation) smaller than Iks- Such rel- 
atively accurate data are necessary for the development 
of approximate methods that can be employed for larger 
or more complicated systems^ They also constitute a 
realistic testbed for present and future path integral tech- 
niques. For comparison purposes, the best known data 
computed by the double virial estimator have a statistical 
error of about 10k b in the low temperature region^ 
The total potential energy of the Nei3 cluster is given 

by 



V tot = J2v LJ (r ij ) + J2v c (r i ), 



(27) 



i<j 



where VLj{rij) is the pair interaction of Lennard-Jones 
potential 



V LJ {n 



V '13 



and V c (i"i) is the confining potential 



V c (ri) = e LJ 



\n - R c 

R, 



211 



(28) 



(29) 



The values of the Lennard-Jones parameters o~lj and clj 
used are 2.749 A and 35.6 K, respectively^ The mass of 
the Ne atom was set to to = 20.0, the rounded atomic 
mass of the most abundant isotope. R cm is the coordi- 
nate of the center of mass of the cluster and is given by 



(30) 



Finally, R c = 2a lj is the confining radius. The role of 
the confining potential V c (rj) is to prevent atoms from 
permanently leaving the cluster since the cluster in vac- 
uum at any finite temperature is metastable with respect 
to evaporation. 
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The optimal choice of the parameter R c for the con- 
straining potential has been discussed in recent work^S If 
R c is taken to be too small, the properties of the system 
become sensitive to its choice, whereas large values of R c 
can result in problems attaining an ergodic simulation. 
To facilitate comparisons, in the current work, R c has 
been chosen to be identical to that used in Ref. [la. 



A. Sampling strategy 

The sampling strategy utilized in the present paper is 
similar to the one employed in Ref. |2(] except for the use 
of parallel tempering^&MiSiSLSLSSi^ to cope with possi- 
ble ergodicity problems. We have utilized a number of 
42 parallel streams, each running a replica of the system 
at a different temperature. For each stream, the basic 
Monte Carlo steps consist in moves of the physical coor- 
dinate X; of an individual particle together with the first 
one quarter of the associated path variables or of the last 
three quarters of the path variables for the respective 
particle. Eq. (27) of Ref. ]24| as specialized for the short- 
time approximation constructed in Section IV. B of the 
same reference, shows that the first one quarter of the 
path variables are associated with Schauder functions, 
whereas the last three quarters are special functions de- 
signed to maximize the asymptotic rate of convergence of 
the path integral method employed. Given the analytical 
differences between the Schauder and the special func- 
tions, one expects that the optimal maximal displace- 
ments for the path variables associated with functions 
from the two categories are different. We mention that 
a poor sampling of the path variables associated with 
the special functions might ruin the quartic asymptotic 
convergence of the path integral method employed. For 
this reason, we attempt to update the path variables as- 
sociated with the Schauder functions and with the spe- 
cial functions separately. The physical coordinate Xj is 
updated together with the Schauder functions. Distinct 
acceptance ratios are computed for the two steps. The 
maximal displacements are adjusted in the equilibration 
phase of the computation so that each of the acceptance 
ratios eventually lies between 40% and 60%. 

The basic computational unit is the pass, defined as 
the minimal set of Monte Carlo attempts over all vari- 
ables in the system. Thus, a pass consists of 2 • 13 = 26 
basic steps. Each Monte Carlo attempt is accepted or re- 
jected according to the Metropolis logici 40 i 41 One defines 
a block as a computational unit made up of forty thou- 
sand passes. The length of the blocks is large enough to 
ensure independence between the block-averages of var- 
ious quantities computed. This independence has been 
checked with statistical tests, as described in Ref.j2CJ As 
opposed to the computation performed in Ref. |2JJ, the 
correlation between block averages of different streams 
has not been tested for independence. The explana- 
tion is that these block averages are correlated by the 
parallel tempering algorithm. However, we have tried 



to minimize this correlation by employing separate ran- 
dom number generators for each streams. These random 
number generators are obtained with the help of the Dy- 
namic Creator packag o 42 i 43 and produce highly indepen- 
dent streams of random numbers, as demonstrated by 
the statistical tests performed in Ref. I20L 

A swap between streams of neighboring temperatures 
has been attempted every 25 passes and it has been ac- 
cepted or rejected according to the parallel tempering 
logic i^Mi^iSSiSLSSi^ Any given stream attempts a swap 
with the neighboring streams of lower and higher tem- 
peratures in succession. Because of this swapping strat- 
egy, the streams of minimum and maximum tempera- 
tures are involved in swaps every 50 steps, only. The 
interval [4, 14] has been divided in 40 equal subintervals 
demarked by 41 intermediate temperatures. Thus, the 
lowest temperature stream has run at a temperature of 
4 - (14 - 4)/40 = 3.75 K. The efficiency of the parallel 
tempering algorithm depends strongly on how much the 
distributions for neighboring temperatures overlap. In 
classical simulations, the width of the overlap is inversely 
proportional to the difference between inverse neighbor- 
ing temperatures. It appears then that the optimal divi- 
sion of the interval [4, 14] involves equally spaced inverse 
temperature subintervals. While not the optimal one, our 
choice of equal temperature subintervals has the advan- 
tage that it provides a smoother heat capacity curve. We 
have monitored the acceptance ratios for all 42 streams 
and found values larger than 60% for all simulations per- 
formed. Thus, the overlap between neighboring temper- 
atures is more than adequate. 

As previously mentioned, besides the acceptance ratios 
of swaps, we have also monitored individual acceptance 
ratios for the Metropolis sampling at each temperature. 
We have ensured that these acceptance ratios are be- 
tween 40% and 60% by automatically adjusting the val- 
ues of the maximal displacements for the path variables 
in the equilibration phase of the computation. Numeri- 
cal experimentation has showed that in order to achieve 
a statistical error of less than lfce for heat capacities, it 
suffices to set the length of the data accumulation phase 
to 100 blocks, for a total of 4 million passes per tempera- 
ture. The equilibration phase has consisted of 20 blocks. 
We have therefore employed a number of data accumu- 
lation passes per temperature equal to the one utilized 
by Neirotti, Freeman, and Doll. This facilitates a direct 
comparison between the two heat capacities estimators 
introduced in the preceding section and the double virial 
estimator. 

The discretization step eo entering the finite difference 
schemes has been set to eo = 2~ 12 . We mention that 
computer experimentation has shown that the numerical 
accuracy of the finite difference schemes is at least one 
thousand times smaller than the statistical error for all 
simulations performed and for a large range of eo . Good 
values for eo are any inverse powers of two between 2~ 18 
and 2~ 8 . 

We conclude this section by commenting on the evalu- 



ation of the errors involved in the determination of heat 
capacities. As opposed to energy estimators, heat capac- 
ity estimators are biased. This is apparent from Eq. J5J). 
In a general setting, let us denote by Xi and Yi the block 
averages of two quantities X and Y and let us define 

^ n 1 n 

X n = -yx t and Y n = - V Y t . 



Given a continuously differentiable function f(x,y), we 
have 

f(X n ,Y n )->f((X),(Y)) 

almost surely, but unless f(x,y) is linear in its variables, 
f(X n ,Y n ) is a biased estimator of f((X) , (Y)). In the 
limit that the variables X n and Y n have small fluctua- 
tions around their expected values, the following approx- 
imation holds: 



100 



[/(x„,y„)-/((x),(r»] : 



—f((X),(Y))(X ri 



(x)) 



d 



+ -f((X),(Y))(Y n -(Y)) 



As such, the mean square deviation for the quantity of 
interest is given by the variance of the quantity 

^-f {{X ),(Y))X n + ^-f((X),(Y))Y n , 

variance that can be evaluated with the (again biased) 
estimator 



1 ™ 



n(n - 1) f-f 



d_ 

dx 



f(X n ,Y n )(Xi — X n ) 



d 



-i 2 



+ —f(X n ,Y n ){Yi-Y n ) 



(31) 



The error bars reported in the present work represent 
twice the square root of the above expression. For the 
heat capacity problem, f(x,y) = x — y 2 and the quan- 
tities Xi and Yi represent block averages of the sec- 
ond order and the first order derivatives of the function 
i? n (x, a;/3, e) around the point e = 1 [see Eqs. ©, (JT^J, 
and HU]. 



B. Numerical results 

A graph of the heat capacity computed with the TT- 
method estimator as a function of temperature is found 
in Fig. ^ for each number of path variables employed. 
The sole exception is an additional run performed with 
a number of N = 127 path variables, which produces re- 
sults virtually indistinguishable (i.e., the differences are 




FIG. 1: Heat capacities (in units of fcs) computed with the 
TT-method estimator as a function of T (in Kelvin) for several 
values of N. The error bars (two times the standard devia- 
tion) are comparable to the thickness of the drawing lines and 
are not plotted. 



smaller than the error bars) from the N = 63 results. 
Therefore, the remainder of the simulations have been 
performed using TV = 63 path variables. Table [I] of Ap- 
pendix A contains the values obtained in the N = 127 
simulation for T = 4, 5, . . . , 14 as well as the associated 
error bars. We believe such values are useful both in the 
design of approximate quantum methods and as a numer- 
ical test for present and future path integral methods. 

The modified TT-method estimator produces results 
similar to the direct TT-method estimator. As shown 
in Fig. |21 the curves for the two estimators are virtu- 
ally indistinguishable. Fig. [21 also contains the classical 
heat capacity as a function of temperature. As appar- 
ent from Fig. the modified TT-method estimator has 
a larger variance in the low temperature region than the 
TT-method estimator. Though they seem to diverge to 
infinity as T — > 0, the error bars of both quantum es- 
timators are comparable to the error bars of the classi- 
cal estimator for the range of temperatures investigated. 
In the low temperature range, the error bars are about 
10 times smaller than those reported by Neirotti, Free- 
man, and Doll 15 for the double virial estimator. Taking 
into consideration that the same number of Monte Carlo 
points has been employed, the TT-method estimator is 
over 100 times more efficient than the double virial esti- 
mator. We mention that the improvement has little to 
do with the path integral technique that has been uti- 
lized. Provided that enough path variables are consid- 
ered, the variance of the estimators is independent of the 
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FIG. 2: Classical heat capacities (cTT) and quantum heat 
capacities by the TT-method estimator (qTT) and the mod- 
ified TT-method estimator (qmTT) as functions of temper- 
ature. On this scale, the curves for the last two quantities 
overlap. The heat capacities are given in units of ks, whereas 
the temperature is given in Kelvin. The number of path vari- 
ables employed for the quantum results is N = 63. The error 
bars (two times the standard deviation) are comparable to 
the thickness of the drawing lines and are not plotted. 




FIG. 3: Error bars (in units of ks) for classical heat capacities 
(cTT) and quantum heat capacities by the TT-method esti- 
mator (qTT) and the modified TT-method estimator (qmTT) 
as functions of T (given in Kelvin) . Also plotted (solid line) is 
the absolute value of the difference Att = (CV) J T — {CV)™ TT 
between the heat capacity values computed with the help of 
the TT-method and modified TT-method estimators. 



path integral technique. At least in one other instance, 
such a significant improvement in the efficiency of a path 
integral technique has been eventually attributed to a 
superior estimatori^SiAi 

As emphasized in Ref. [53, the agreement between the 
T-method and the H-method energy estimators consti- 
tutes an important test for the convergence of the path 
integral methods. The heat capacity analog is repre- 
sented by the agreement between the TT-method and 
the TH-method estimators. The latter estimator is ob- 
tained by temperature differentiation of the H-method 
energy estimator. The temperature differentiation can 
be performed by finite difference in a way similar to 
the present implementation of the TT-method estima- 
tor. However, the evaluation of the H-method estimator 
requires knowledge of the first order derivatives of the 
potential. Since these derivatives have been computed 
anyway in the modified TT-method estimator simulation, 
we have also evaluated the H-method energy estimator in 
the respective simulation. A temperature differentiation 
with the help of the formula 



(a 



TH 



(E) 



(E) 



Pi+i - A- 



(32) 



has produced the curve in Fig. 0] figure that also plots 
the TT-method heat capacity estimator, for comparison. 
The agreement between the two curves is surprisingly 
good. In fact, the maximum difference between the two 
curves is about 1.5/c.b, a value that is comparable to the 
error bars achieved in the present simulations. 

We say that the agreement is surprisingly good be- 
cause several factors concur against such an agreement. 
First, numerical differentiation of Monte Carlo data is in 
general a difficult task, unless the data at different tem- 
peratures are strongly correlated so that the resulting 
curve is smooth. In this regard, the parallel tempering 
technique is of great help because it brings the necessary 
correlation into the simulation. From the quality of the 
numerical differentiation, we estimate that the correla- 
tion is substantial. For instance, if the runs at different 
arbitrarily close temperatures are independent, the re- 
sulting curves fail to be continuous. If the correlation 
is of the type appearing in a Brownian motion, the re- 
sulting curves are continuous but not differentiable. In 
order for the curves to be differentiable, the correlation 
must be even stronger. Though such a strong correla- 
tion has been previously reported^ we are not aware 
of any mathematical or numerical analysis attempting to 
quantify the strength of the parallel tempering correla- 
tion between averages at different temperatures. In the 
light of the application just presented, we believe such an 
analysis would be well justified. 

Second, the temperature step in the numerical differen- 
tiation is significantly larger than the step we have em- 
ployed for the TT-method estimator. Fortunately, the 
quantum effects are strong and the dependence of the 
ensemble energy with the temperature is smooth. As a 
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FIG. 4: Quantum heat capacities in units of ks by the TT- 
method estimator (qTT) and by the TH-method estimator 
(qTH), respectively. The temperature is given in Kelvin. On 
this scale, the two curves overlap almost perfectly. The max- 
imum difference between corresponding values on the curves 
is about 1.5fcs. 



consequence, the accuracy of the finite-difference scheme 
is comparable to the statistical errors. 

A third factor that could prevent an agreement be- 
tween the TT-method and the TH-method heat capacity 
estimators is the lack of convergence of the path integral 
method employed. The agreement provides additional 
evidence that N — 63 path variables are sufficient for the 
range of temperature studied and for the path integral 
technique utilized. 

Yet a fourth reason for disagreement is poor Monte 
Carlo sampling. Energy estimators are unbiased estima- 
tors, as opposed to heat capacity estimators, which are 
biased. As a consequence, energy estimators and heat 
capacity estimators generally have different sensitivities 
to the quality of the sampling, with the latter ones be- 
ing more sensitive to quasiergodicity problems. This may 
result in disagreement between the heat capacities com- 
puted with the help of estimators and the ones computed 
by using energy differences of the type given by Eq. (|32H . 

We conclude this section by noticing that the high- 
temperature part of the quantum heat capacity plotted 
in Fig. 121 does not coincide with the results reported in 
Ref. [la. The cause of this difference is the fact that 
Neirotti, Freeman, and Doll have mistakenly utilized a 
confining potential with a radius R c = Actlj instead of 
2<7£j, the value they have reported. 



IV. SUMMARY AND CONCLUSIONS 

The main result of the present paper is the finding 
that the evaluation of the main thermodynamic prop- 
erties of a quantum canonical system, namely average 
energy and heat capacity, can be performed in a fast and 
reliable fashion without calls to first or second deriva- 
tives of the potential. This can be accomplished by a 
finite-difference scheme applied to the T-method energy 
estimator and TT-method heat capacity estimator, re- 
spectively. The derivation of these estimators is rather 
trivial, consisting of simple temperature differentiations 
of the partition function. As emphasized in the intro- 
duction, the key observation is that the Feynman-Kac 
formula and its finite-dimensional approximations must 
be written in a form with the temperature dependence 
of the paths buried into the potential. Such a transfor- 
mation is possible for all path integral techniques and it 
should constitute the starting point for the derivation of 
various energy and heat capacity estimators. 

We have also proposed an analytical heat capacity es- 
timator, called the modified T-method estimator, that 
might prove useful whenever the first derivatives of the 
potential are available. However, this estimator has a 
slightly worse behavior at low temperature than the di- 
rect TT-method estimator and may be quite expensive 
for large-dimensional systems because of the quadratic 
scaling of the number of path integrals that must be com- 
puted with the dimensionality of the system. For exam- 
ple, in the case of the Nei3 cluster, the code based on the 
modified heat capacity estimator has been 50% slower 
than the code utilizing the finite-difference scheme. 

The heat capacity estimators utilized in the present 
paper have favorable variances when compared to the 
double virial estimators. This has been clearly demon- 
strated for a Lennard- Jones realization of Nei3, a realis- 
tic physical system that is representative of many other 
applications. To the authors' knowledge, the heat capac- 
ities results obtained for the Nei3 cluster are the most 
accurate to date. 
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APPENDIX A: TABLE OF HEAT CAPACITIES 
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TABLE I: Heat capacities and error bars of the Nei3 cluster as functions of temperature. A number of N = 127 path variables 
have been utilized. The error bars are two times the standard deviation. The temperature is measured in Kelvin, whereas 
the heat capacities are given in units of ks- The heat capacity pick value, as obtained by maximizing a cubic spline function 
interpolating the computed data, is 74.47 ± 0.54 ks and is attained at the temperature of T pca k = 8.97 K. 
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5 
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9 


(Cv) 


8.26 ± 0.80 


11.81 ± 0.52 


18.14 ± 0.44 


29.30± 0.56 


54.74 ± 0.86 


74.45 ± 0.54 


T 


10 


11 


12 


13 


14 




(Cv) 


61.70± 0.51 


48.87± 0.36 


43.05 ± 0.30 


40.78 ± 0.23 


40.09± 0.27 
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